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Abstract 

Three  related  successive  approxima¬ 
tion  schemes  for  determination  of  optimal 
trajectories  are  developed  with  particular 
attention  to  treatment  of  inequality  con¬ 
straints  on  control  variables.  The  first 
is  a  gradient  method  based  upon  a  Eucli¬ 
dean  metric  with  appropriate  modification 
for  handling  of  inequalities;  the  second 
employs  a  Min  operation  without  use  of  a 
metric;  and  the  third  features  a  special 
integral -absolute  value  metric.  The  Pon- 
tryagin  Principle  is  employed  for  con¬ 
struction  of  successive  control  function 
approximations.  All  three  schemes  employ 
an  adjoint  system  for  computation  of  in¬ 
fluence  functions  and  a  'penalty  function' 
technique  for  handling  of  constraints  on 
terminal  values. 

Illustrative  calculations  are  pre¬ 
sented  for  planar  Earth -Mars  transfer 
with  rocket  thrust  variable  between  lim¬ 
its  T^^  <  T  <  T2.  The  relative  merits  of 

the  techniques  are  discussed  from  the 
viewpoint  of  digital  computation. 

Introduction 

Recent  work  on  successive  approxima¬ 
tion  techniques  for  numerical  solution  of 
variational  problems  involving  differen¬ 
tial  equations  as  subsidiary  conditions 
has  provided  a  clear  indication  of  the 
practical  usefulness  of  this  class  of 
methods  when  employed  in  conjunction  with 

high  speed  digital  computatlonl> 2,3,4 , 

the  present  paper  we  report  some  recent 
developments  of  techniques  which  are  ap¬ 
plicable  to  problems  featuring  inequality 
constraints  imposed  on  the  control  varia¬ 
bles  . 


Problem  Formulation 

The  basic  system  of  differential 
equations  describing  the  process  of  inter¬ 
est  is  presumed  given  in  first  order  form 


Xf  -  x^,  y^,  -  -,  y^,  t)  (1) 

i  -  1,  -  -  -,  n 


The  variables  Xj^,  i  -  1,  -  -,  n,  are 
state  variables  and  the  variables 
yk,  k  *  1,  -  are  control  variables. 

The  functions  gj^  are  assumed  to  possess 
continuous  first  partial  derivatives  with 
respect  to  their  arguments. 

Initial  conditions  x^(tQ)  ••  x^^, 

i  -  1,  -  -,  n,  are  prescribed.  In  addi¬ 
tion  some  of  the  terminal  values  Xj^(tf)  ■ 

x<  may  also  be  prescribed.  The  terminal 
f 

time  t£  will  usually  not  be  fixed;  thus 
its  value  is  regarded  as  'open'  for  opti¬ 
mization  purposes. 

We  wish  to  minimize  some  function 

P  »  P(x,  ,  -  -,  X  ,  t.)  (2) 

if  n^  r 


of  the  terminal  values  of  the  state  varia 
bles  and  the  terminal  time.  This  is  the 
general  problem  of  Mayer,  encompassing  a 
large  class  of  current  problems  in  flight 
performance,  control,  and  guidance.  It 
has  recently  been  observed  by  Hoelker^ 
that  similarities  in  mathematical  treat¬ 
ment  of  these  problems  provide  a  strong 
influence  toward  unification  of  future 
efforts . 


'Penalty  Function'  Treatment  Of 
Terminal  Constraints 

We  may  relax  the  requirement  for 
meeting  fixed  terminal  conditions  on  the 
X£  in  favor  of  an  approximation:  the 
addition  of  terms  to  the  function  P  of 
the  form 
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m 


n 


n 


i 


I 

P 


p  +  i 


(3) 


1-1  k-1 


(6) 


Kj  >  0 

(The  summation  ranges  only  over  those 
variables  subject  to  specified  terminal 
conditions .) 

The  notion  underlying  this  approxima¬ 
tion  is  that  the  solution  of  the  problem 
Min  p'  will,  under  appropriate  condition^ 
tend  toward  the  solution  of  the  original 
problem  as  the  magnitudes  of  the  Kj  are 
increased.  This  idea  is  due  to  Courant^ . 
Its  basis  and  application  are  discussed  in 
some  detail  in  Ref.  3.  We  may  now  deal 
with  the  problem  of  minimizing  P’  with 
terminal  values  of  the  'open.' 

Adjoint  System 

1  O  -J 

As  in  earlier  work^'-*'  ,  we  employ  an 
adjoint  system  of  differential  equations 
in  variables  obtained  from  the  line¬ 

arized  version  of  (1) 


which  may  be  verified  directly  by  differ¬ 
entiation  and  use  of  (4)  and  (5) .  Inte¬ 
grating  both  left  and  right  members  of  (6) 
between  definite  limits  we  have 


t  i-1  k-1 


o 


In  terms  of  a  function  H  defined  by 
H  (x^,  -  -,  x^,  -  -,  X^, 


6x 


i 
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J-1 


I 


k-1 


(4) 


y^,  -  t) 


(8) 


i-1,  -  -,  n 


[The  linearization  is  in  the  neighborhood 
of  yi^  =  3  3  solution  of 

system  (1) .  ] 

To  obtain  the  adjoint  system,  one  discards 
the  control  terms,  transposes  the  matrix 
of  coefficients,  and  changes  the  signs  of 
the  right  members: 


n 


j-1 


this  may  be  written: 


(9) 


The  partial  derivatives  are  func- 

i 

tions  of  t  only,  being  evaluated  along 

yk  ”  yk<t),  -  Xi(t)  . 

The  most  significant  property  of  ad¬ 
joint  systems  is  given  by 


We  shall  need  a  version  of  (9)  which 
incorporates  higher  order  effects,  in  par¬ 
ticular  one  for  which  linearization  with 
respect  to  the  control  variables  is  a- 
voided.  Our  development  from  the  proper¬ 
ties  of  adjoint  systems  follows  that  of 
Ref.  7.  Working  temporarily  with  finite 
increments  Ax^  =  Ayj^  “  yk  '  ^k’ 

we  get  from  (1) 


2 


AXi  -  x^  +  Ax^,  yj^  +  Ay^^,  -  ->  y^  +  Ay^,  t) 


8j^(X]_i  “  “*  x^,  y^,  -  *,  y^,  t) 


(10) 


From  integration  of 


n 


i-1 


^i  ^’'i  ^ 


Ax^ 


i-1 


(11) 


we  then  obtain 


Xj^  Ax^  dt 


n 

>'i[8i(Xl  +  AXj^,  -  -,  +  Ax^, 

i-1 


-  gi(Xi»  -  x^,  y^,  -  -,  y^,  t) ]  dt  (12) 

If  we  employ  an  expansion  of  the 
functions  gj.  in  a  Taylor  series  in  the 
Axj^,  including  a  remainder  term 


gj^(x^  +  Axj^,  - 

“  g^(x2^» 


x^  +  Ax^,  yj^  +  Ayj^,  -  y^  +  Ay^,  t) 

■  •'  *n'  yi  +  ^^1'  ■  y/  +  ^y/' 


n 


j-i 


■'  *n'  yi  ^yi’  ■  ■'  y^  ■*■  ^yf* 


z 

j-l  s-l 


- 

(^1  +  ^n  ^n^’^n' 

j  8 

7l  +  Ay^,  -  -.  Vl  +  Ay/,  t)AxjAx8 


(13) 
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where  0  <  ?q  <  1,  q  -  1,  -  n,  then 
(12)  becomes 
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i-1 


"f  "f  n 


t  1-1 
o 


-  ■  ■'  ’‘n’  ^1'  ■  ■'  '^V  ^ 


.  f  n  n 


t  i-1  j  =  l 
o 


X. 

1 


_  _  _ 

Fx"  ^’‘l’  ■  ■'  ^n'  ^1  +  ^^1'  •  •' 


-  -  - 
ET  ^’^1'  ■  ■'  ’^n'  yi'  •  ■' 


AXj  dt 


+  i 


■f  n  n  n 


t  i-1  i-1  s-1 

o  ■’ 


n2 

S  gj^  _ 

N — =: —  (x,  +  ^,Ax,  ,  -  X  +  ?  Ax  , 
i  dx.ox  '  1  11'  '  n  n  n' 

J  s 


Yl  +  Ay^^,  -  +  Ay^,  t)  Ax.Ax^  dt  (14) 


Note  that  the  series  development  assumes 
the  existence  of  the  second  partial  deriv¬ 
atives  of  the  g£  with  respect  to  the 
state  variables.  In  terms  of  the  function 
H*  defined  by  (8) ,  the  expression  (14) 
may  be  written; 


AXj^ 


^f  „^f 


t  o 

o 


★ 


[H  (yj^  +  Ayj^,  •  y/  +  Ay^,  t)  -  H  (y^^,  -  y^,  t)  ] 


dt 


in  which  only  the  term  of  the  right  member 
dominant  for  small  Ay^  is  shown.  Since 
the  functions  Xj^(t)  ,  ^^(t)  may  be  re¬ 
garded  as  known  functions  of  t  evaluated 
from  (1)  and  (5)  for  the  control  functions 
yi^(t) ,  the  function  H*  appearing  in 
(15)  may  be  considered  a  function  of  the 
Ayit  and  t  only. 


Terminal  Time  Criterion 

We  employ  as  trajectory  termination 
criterion  the  vanishing  of  the  time  de- 


I 

rlvatlve  of  the  function  P  (xj^, 


n' 


t)  - 


n 


1-1 


(15) 


-.Xn.t): 


(16) 


This  amounts  to  a  one -dimensional  search 
for  the  value  of  tf  which  minimizes 

tf) ,  performed  concur¬ 
rently  with  numerical  Integration  of  the 
basic  system  (1) .  The  terminal  time  tf 
is  determined^by  a  zero  of  (16) ,  t  >  t^, 

for  which  p'  is  negative, 

dt^ 


Effect  Of  Control  Variations  On  The 
Function  P' 


Beginning  with  some  first  approxima¬ 
tion  for  the  control  functions  yj^  - 
and  the  corresponding  solution  of  (1) , 

Xf  »  Xf (t) ,  we  wish  to  obtain  successive 
decrements  in  the  performance  function 
P',  eventually  approaching  a  minimum.  We 
assume  that  this  first  solution  of  (1)  is 
terminated  according  to  (16) ,  and  we  des¬ 
ignate  the  terminal  time  tf  -  tf.  If  the 
control  functions  are  altered,  producing 
increments  AXf(t)  in  the  state  variables, 
a  first  order  approximation  to  the  incre¬ 
ment  in  the  terminal  value  of  P'  is 
given  by 


t 

AP 


n 


Ax^(tp 


n 


Xi(tf) 


At^  (17) 


f 

AP 


[H*(yj^  +  Ayj^, 

t 

o 


y/  +  t) 


-  H*(yj^,  -  -,  y^,  t)  1  dt  (19) 


the  higher  order  terms  having  been 
dropped. 


Modified  Gradient  Method 

Sketching  first  the  basic  gradient 
method  for  obtaining  negative  Increments 
in  P',  we  en^loy  (19)  for  construction 
of  control  variations  Ayjj(t)  .  We  wish  to 
limit  the  'step  size'  in  some  way  as  to 
guarantee  the  predominance  of  first  order 
effects.  A  constraint  on  the  'distance' 
covered  in  the  step  as  measured  in  the 
Euclidean  sense 


Ayj^^  dt  -  ,  k  -  1,  -  -,  i  (20) 

*^t 

o 

is  the  conventional  means  employed.  If 
the  constants  ak  are  taken  sufficiently 
small,  then  (19)  will  be  a  good  approxima¬ 
tion  to  the  actual  change  in  P'  and  it 
is  sensible  to  seek  a  minimum  of  (19)  sub¬ 
ject  to  (20)  . 


in  which  the  partial  derivatives  of  P' 

are  evaluated  for  xj  "  xi  ,  tf  •  tf. 

f_  f 

But  by  our  choice  of  tf  determined  by 

(16)  ,  the  second  member  on  the  right  of 

(17)  vanishes. 

Referring  to  (9)  and  (15)  and  noting 
that  the  Axf(tQ)  are  zero  as  a  conse¬ 
quence  of  fixed  initial  conditions  on  the 
Xf,  we  see  that  if  boundary  conditions 

\('£)  ■  ir  <*») 

^  c. 


are  imposed  upon  the  adjoint  variables 
evaluated  along  yk(t),  xf(t),  then  the 
following  form  of  (15)  is  obtained: 


For  this  purpose  we  apply  the  Pontry- 
agin  theory7,8,9^  the  full  advantage  of 
this  approach  becoming  clear  as  we  later 
proceed  with  treatment  of  inequality  con¬ 
straints  on  the  y^.  The  integrals  (19) 
and  (20)  are  represented  as  the  terminal 
values  of  variables  zf(t),  -  -,  z/+i(t) 
satisfying  the  system: 


^k  "  ^yk^  '  ^k^‘^o^  “  ®  “  ®k^ 

k  -  1,  -  -  -,  /  (21) 

^£+1 "  »*<yi  ^yi'  •  y£  +  *=> 

•ff  _ 

(y^,  ~  ",  y^,  t)  , 


5 


(22) 


Thus  a  statement  in  terms  of  an  auxiliary 
Mayer  problem  is  obtained:  Min  z^_|  j^(t£)  . 

Introducing  adjoint  variables  Pi(t), 
■  Pi+i(t)  and  a  function  H  defined 
by 

H  -  ^  Pk  Ayk^ 

k-1 

+  Pi+i[H  (y^  +  Ayj,  ■  Yi  +  Ay^,  t) 

•  H*(yj^,  -  y^,  t)  ]  (23) 

we  then  write  the  conditions 

Pk  -  0  ,  k  -  1,  -  i  (24) 

Pm  -  0  ,  P^+i(t£)  -  1  (25) 


H(Ay£,  -  Ay^,  Pj^,  -  p^^^,  t) 

A  ^ 

^  H(Ay£>  *  Ay^>  Pl^  *  P/+l^  (2^) 


A  _ 

The  Ayk  “  Ayk(t)  minimizing  zi+i(tf) 
are  determined  from  (26) ,  the  Pontryagin 
principle,  whose  alternate  form  is 

Min  H  (27) 

Ay^,  -  Ay^ 

The  expressions  (26)  and  (27)  are 
equivalent  to  the  necessary  condition  of 
Welerstrass,  and  represent  a  generaliza¬ 
tion  thereof  if  inequality  constraints  on 
the  Ayk  are  operative,  as  will  later  be 
the  case7>9.  in  the  absence  of  such  con¬ 
straints,  the  minimum  of  H  sought  in 
(27)  must  also  be  a  stationary  point, 
this  following  from  the  differentiability 
property  of  the  functions  gi  assumed 
earlier.  The  vanishing  of  the  partial 
derivatives 


then  determines  the  control  increments  as 


1  5H* 

2Pk  ^^k 


k  -  1,  -  -,  i  (29) 


The  constants  pi,  -  -,  p^  must  sat 
isfy  the  requirements 


2k(tf) 


(30) 


or 


Pk  ' 


k  -  1,  -  t  (31) 


Evidently  the  negative  sign  corresponds  to 
a  maximum  of  zi+i(tf)  •  AP'  and  the  pos¬ 
itive  sign  to  a  minimum.  Choosing  the 
latter,  we  obtain 

i, 

2 


(32) 


We  have  found  the  direction  of 
'steepest  descent'  in  the  sense  of  the 
Euclidean  metric  (20) ,  and,  according  to 
(32),  the  increment  AP'  must  be  negative 

(t)  evalua- 

Syk 


or  zero.  The  functions 


ted  for  yk  “  yk(b)  assume  the  role  of 
gradient  components  and  the  ' steepest 
descent'  process  is  governed  by 


i,  -  I 


^yk 


bk<0 


(28) 
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k  -  1, 


£ 


(33) 


We  now  turn  attention  to  means  for 
handling  Inequality  constraints  on  the 
control  functions  y)(.  It  Is  assumed 
that  the  Inequalities  can  be  converted  to 
the  form 


^ki  <  yfc  <  yk2  ' 
k  -  1.  -  i  (34) 


which  Is  usually  the  case  In  applications. 
The  development  up  through  (26,  27)  Is 
Identical  for  this  case,  except  that  the 
Min  operation  In  (27)  oust  be  performed 
subject  to  the  additional  constraints 
(34); 


Min  H  (35) 

Ay^,  •  Ay^ 


^ki  <  ^k  <  yk2 


An  examination  of  the  minimum  prob¬ 
lem  (35)  for  vanishingly  small  step  size 
|aj(|  leads  to  the  formulas 


<  0  (36) 


If  y](  Is  an  Interior  point  of  the  In 
terval,  l.e.,  for 


yki  <  ^k  <  yk2 


At  the  lower  limit  yjj  ■■  yj^^  , 


A  w 

\  ^ 


Ay.  -  0 


If 


5^<o' 

5yv  - 


\b^  <  0  (36) 


while  at  the  upper  limit  yj^  ■  yy^2 


Ay,  -  0 


”  '’k  I?;; 


If 


If 


^  *  N 

^<0 

Sy^  A 


W  <  0  (39) 


^  * 

|S->  0 

*>'k-  ! 


Thus  we  arrive  at  a  modified  gradient 
method. 


In  a  descent  process  for  which  the 
steps  are  vanishingly  small,  essentially 
a  continuous  process,  the  governing  equa¬ 
tions  (36) ,  (38) ,  (39)  will  succeed  In 
holding  the  control  variables  yj^  within 
the  desired  range  yj^j^  <  yk  yk2’  There 
will  be  difficulty,  however,  if  finite 
steps  are  taken,  and  to  avoid  this,  It 
will  be  necessary  to  'tiim'  the  functions 
yjj  +  Ayjj  obtained  for  finite  bj^  to  con¬ 
form  to  the  inequality  (34)  before  employ¬ 
ing  them  in  numerical  integration  of  the 
basic  system  (1) . 


The  general  mode  of  operation  for 
gradient  calculations  has  been  described 
In  earlier  publlcatlonsli 3  for  cases  In 
which  only  a  single  control  variable  is 
employed.  It  consists  of  numerical  solu¬ 
tion  of  the  basic  system  (1)  for  a  number 
of  values  of  the  descent  parameter  b, 
the  parameter  being  varied  systematically 
In  one -dimensional  search  for  a  minimum 
of  p' .  In  problems  Involving  a  number 
of  control  variables,  an  Identical  treat¬ 
ment  may  be  given  to  each  in  turn;  or  al¬ 
ternatively,  a  more  complicated  multi¬ 
dimensional  search  versus  the  bj^  may  be 
Implemented. 


A  Successive  Approximation  Scheme 
Employing  the  Min  Operation 

A  shortcoming  of  a  gradient  process 
is  that  over  intervals  in  which  is 

^yk 

small  In  magnitude,  the  corresponding 
changes  In  y;^  will  be  small.  After  sev¬ 
eral  steps  the  yj^  may  still  be  far  from 
their  optimal  values  over  such  'insensi¬ 
tive'  intervals  owing  to  this  feature  of 
the  gradient  process.  This,  of  course, 
stems  from  the  rather  arbitrary  Imposition 
of  the  Euclidean  distance  measure.  From 
an  engineering  viewpoint  this  Is  unimpor¬ 
tant  If  only  the  value  of  P  is  of  inter- 
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est,  as  Is  often  the  case  in  flight  per¬ 
formance  work.  It  is  inconvenient  if  a 
family  of  neighboring  extremals  are  re¬ 
quired,  as,  for  example,  in  connection 
with  a  guidance  study,  because  this  re¬ 
quires  that  additional  computations  be 
performed  to  converge  the  control  variable 
histories  to  within  the  desired  accuracy. 

It  has  been  speculated  that  an  ap¬ 
propriate  procedure  for  treatment  of  this 
situation  is  transition  to  a  scheme  for 
systematic  numerical  solution  of  the 
Euler  equations,  and  this  appears  plausi¬ 
ble  on  first  consideration.  One  finds, 
however,  that  the  appropriate  linear  com¬ 
binations  of  adjoint  solutions  do  not 
yield  a  good  approximation  to  the  multi¬ 
plier  functions  of  the  'indirect'  theory, 
and,  in  particular,  the  required  initial 
values  of  the  multipliers  may  be  suffi¬ 
ciently  in  error  to  cause  difficulty  in 
an  iterative  adjustment  process. 

The  successive  approximation  scheme 
of  the  present  section  was  originally  de¬ 
veloped  for  refinement  of  control  pro¬ 
grams  of  near-minimal  gradient  solutions 
into  close  approximations  to  Euler  solu¬ 
tions,  i.e.,  to  accelerate  convergence  in 
the  later  stages  of  the  descent  process. 

It  appears,  however,  to  possess  merit  as 
a  primary  computational  scheme,  as  will 
be  discussed  later  in  connection  with  an 
example . 

Instead  of  employing  the  somewhat 
arbitrary  Euclidean  metric  (20) ,  we 
choose  an  equally  arbitrary  alternative: 
discard  the  use  of  a  metric,  calculating 
new  control  variables  yj^*  from  the  op¬ 
eration 

Min  H*(y)  (40) 

^kl  <  ^k  <  yk2 

k  -  1,  I 

In  adopting  the  control  yj^  -  yk*(0 
generated  by  the  Min  operation  as  our  next 
approximation,  we  risk  the  violation  of 
our  linearizing  assumptions,  for  this  may 
represent  a  large  step  process.  More  con¬ 
servatively,  we  may  elect  to  replace  the 
large  step  by  an  exploratory  series  of 
small  ones,  setting 

Xk  “  ^  ^ 

k  =  1,  -  -,  i 


and  evaluating  p'  versus  a  one¬ 

dimensional  search  analogous  to  that  per¬ 
formed  versus  step  size  (the  bk)  in 
gradient  computations . 

There  is,  of  course,  a  question  of 
convergence  with  this  scheme,  hinging  on 
whether  or  not  the  increment  AP'  given 
by  (19)  is  negative  for  some  C  in  the 
range  0  <  C  <  1  normally  explored.  Evi¬ 
dently  a  sufficient  condition  for  AP'  <  0 
for  small  ^  is 


iim[H*(yj^+C(yi^*-yp,  -  -,  y * -y i) ,  t) 

?->o 


-  H  (y^,  -  -,  y^,  t)  ] 


SH*.- 


k=l 


X^.  t) 


*  - 


(Xk  -  Xk)  < 


(42) 


which  requires  that  the  directional  de¬ 
rivative  of  H*  be  negative  in  the  direc¬ 
tion  of  the  minimum  point  yk  “  Xk*'  This 
requirement  will  be  met  globally,  i.e.  for 
all  admissible  starting  points  yk  =  Xk* 
if  the  function  H*(yi,  -  -,  y/,  t)  pos¬ 
sesses  no  stationary  points  or  interior 
extrema  in  the  region  defined  by  (34) , 
other  than  the  minimum  at  yk  =  Xk'’*^' 

In  the  original  conception  of  this 
technique  as  a  refinement  scheme,  only  a 
local  version  of  the  requirement  discussed 
in  the  preceding  paragraph  requires  con¬ 
sideration,  and  this  requirement  is  auto¬ 
matically  met  if  the  gradient  process  has 
progressed  sufficiently  before  transition 
to  the  Min  H*  scheme.  From  the  viewpoint 
of  more  general  applications,  the  require¬ 
ments  on  the  function  H*  discussed  above 
constitute  a  limitation;  yet  the  class  of 
problems  for  which  convergence  is  assured 
a  priori  is  sufficiently  large  as  to  war¬ 
rant  considerable  interest  in  the  scheme. 


A  Successive  Approximation  Scheme 
Employing  An  Integral -Absolute 
Value  Metric 

A  third  method  for  computing  optimum 
trajectories  by  successive  approximations 
is  obtained  by  substituting  an  integral- 
absolute  value  metric  for  the  integral 
square  metric  of  (20) : 


(48) 


i 

/ 

lAykl  dt  -  >  0 

*^0 


Pk  -  0  ,  k  -  1,  -  i 

Pi+1  “  ^  ^ 


k  -  1,  -  i  (43) 


Although  apparently  not  limited  to  cases 
in  which  the  optimal  control  is  bang-bang, 
this  method  appears  to  offer  the  advantage 
in  such  cases  that  if  the  optimal  control 
is  bang-bang  and  the  first  approximation 
is  taken  to  be  bang-bang,  all  successive 
approximations  will  also  be  bang-bang. 

By  applying  the  Pontryagln  theo- 
in  a  manner  similar  to  that  of 
(21)  to  (27) ,  we  obtain: 

Zk(tf)  "  k  -  1,  -  i  (44) 

and 

*i+l  “  H*(yj^+Ay^,  -  -,  y^+Ay^,  t) 

H  (y^^j  “  y t)  ,  ^_g+i(^Q)  “  ®  (45) 


From  (19)  the  function  we  wish  to  minimize 
is 

AP'  -  2i+l(tf)  (^6) 


The  Ayk  are  now  considered  as  new  con¬ 
trol  variables  with  time-varying  con¬ 
straints  which  are  compatible  with  the 
constraints  on  the  control  variables  yk- 
As  in  the  previous  cases,  the  constraints 
on  yk  are  assumed  to  be  of  the  form 

^ki  <  ^k  <  yk2 

A  necessary  condition  for  AP '  to  be 
a  minimum  is  that  the  Ay^  minimize  H 
for  all  to  <  t  <  t£  subject  to  (49) . 

This  is  expressed  as 


Min  H 

Ay^,  -  Ay^  (50) 


Convergence  is  assured  by  limiting  the 
magnitude  of  the  ak  defined  in  (43) . 
However,  we  will  see  that  this  no  longer 
requires  that  the  change  in  the  control 
variable  be  small  at  any  particular  t. 
From  (48)  and  (49)  we  obtain 

i 

H  “  Pk^'^^k^  ”*(yi+Ay£.  - 

k=l 

y^+Ay^,  t)  -  H*(yj^,  -  -,  y^,  t)  (51) 


Thus  we  have  a  new  mininum  problem  to 
which  we  may  apply  the^Pontryagln  princi¬ 
ple.  A  new  function  H  is  defined  by 

t 

Pk'^yki  Pmf«*<yi+^yi'  - 

k-l 

■  *^*(^1'  ■  ■'  y^'  t)  ]  (47) 


where  Pk  is  a  constant  (pk  =  Pk)  ■ 

Since  the  ak  were  not  previously  speci¬ 
fied,  aside  from  the  requirements  that 
they  be  positive  and  small  in  magnitude, 
we  may  consider  the  Pk  as  constants  to 
be  later  found  by  performing  a  search, 
evaluating  AP'  versus  pk- 

A 

The  function  H  to  be  minimized 
given  by  (52)  can  be  thought  of  geometri¬ 
cally  as  the  superposition  of  a  cone  on 
the  surface  AH”: 


where  Pi(t) ,  -  -,  pi4.i(t)  are  adjoint 
variables  satisfying  the  system  of  differ¬ 
ential  equations 


AH  =  H  (y^+Ay^, 


-  H''(yi,  -  -,  yn,  t) 


(52) 
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Let  us  now  consider  the  treatment  of  one 
control  variable  at  a  time .  We  rewrite 
(52)  as 

H  -  H*(yj+Ayj,  -  y^,  t) 

-  H*(yj,  -  y^,  t)  (53) 

Examining  (53) ,  we  see  that  the  minimum  of 
H  with  respect  to  Ayj  occurs  either  at 

a  stationary  point  (SH/SAyj^  =  0)  ,  at  the 
apex  of  the  cone  where  Ayj^  =  0,  or  on 
the  boundaries  determined  by  (49)^,  For 
pj  very  large,  the  minimum  of  H  occurs 
at  the  apex  of  the  cone  for  all  time  and 
Ayk  =  0  everywhere .  As  pj  is  made 
smaller,  a  point  is  reached  where  the  min¬ 
imum  changes  from  the  apex  of  the  cone  to 
a  boundary  or  an  interior  point  over  an 
infinitesimal  interval  of  time.  This 
value  of  pj  we  will  designate  as  the 
threshold  value  of  pj .  If  we  continue  to 
decrease  pj  a  small  amount  below  its 
threshold  value,  Ayk  will  take  on  non¬ 
zero  values  over  finite  intervals  of  time. 
It  is  to  be  noted  that  the  Ayk  at  any 
time  are  themselves  not  necessarily  small. 
The  final  selection  of  pj  is  determined 
by  evaluating  P'  via  (1)  for  several 
values  of  pj  below  its  threshold  value 
and  fitting  a  parabola  to  the  points.  The 
minimum  of  this  parabola  then  determines 
the  optimal  Pj . 

To  illustrate  the  possible  advantages 
of  employing  the  integral -absolute  value 
metric,  let  us  consider  a  system  linear  in 
the  yk  and  for  which  the  optimal  control 
is  bang -bang.  For  this  case  H*  can  be 
written  in  the  form 

H*  “  ‘t>j^(x.  A,  t)  +  yk‘f’2(’'»  t)  (54) 

A 

giving  for  H: 

H  -  Pkl^yk^ 

If  the  first  approximation  to  the  control 
function  is  bang-bang,  it  is  easily  seen 
that  all  successive  approximations  are 
also  bang-bang  since  this  property  propa¬ 
gates,  due  to  (55) . 


tained  when  Pk  is  equal  to  the  maximum 
value  of  |'t>2(t)  |  in  a  region  where  it  is 
possible  to  change  yk  subject  to  (55) 
and  (49) .  If  we  decrease  Pk  a  small 
amount  beyond  this  threshold  value,  Ayk 
will  change  by  its  maximum  value  for  a 
limited  period  of  time. 


Low  Thrust  Example 

For  illustrative  computations  we 
choose  the  problem  of  planar  low  thrust 
transfer  between  circular  orbits,  which 
has  been  employed  in  earlier  technique  de- 
velopments3, 10 ,  xhe  system  of  equations 
governing  the  motion  is  given  by: 


Radial  Acceleration 


V  ■  A 

R  o 


-r-  +  -  sin  0  (56) 


Circumferential  Acceleration 


uv  ,  T  a 
-  -S“  +  —  cos  0 
R  m 


(57) 


Radial  Velocity 

R  =  =  u  (58) 


Circumferential  Angular  Velocity 


(59) 


Propellant  Expenditure 

m  -  g3  “  -  ^  (60) 

e 


The  exhaust  velocity  Vg  is  taken 
constant  and  the  thrust  T  variable  be¬ 
tween  fixed  limits: 


Tj^  <  T  <  T2  (61) 

If  we  assume  a  fully  throttlable  rocket, 
with  a  control  parameter  q. 


0  <  q  <  1 


For  Pk  very  large  (compared  to 
|')>2(t)l)>  '^yk  is  equal  to  zero  for  all 
t.  The  threshold  value  of  Pk  is  ob- 
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(62) 


then 


(68) 


T  -  T2il 


(63) 


The  function  H  takes  the  form 


*  /  2  IR 

»  ■  MT 


T2TI 

+  -  sin  6 

ra 


■  -»  I  uv  ,  2  ‘  Q 

+  “  "S"  +  -  cos  9 

2  \  R  m 


ToTI 

+  A3U  +  ^  -  X5  — 

e 


(64) 


and  the  adjoint  variables  satisfy  the 
system  of  differential  equations 


i  -  1,  -  5 

[Note  that  for  an  optimal  trajectory  the 
function  H*  becomes  Che  Hamiltonian 
function,  i.e.  as  the  successive  approxi¬ 
mation  process  converges,  H*  -♦  H.  ] 

The  function  p'  is 


p’  -  P  +  MKj^(Uf-u)^  +  K2(v£-v)2 

^2  ^2 

+  K3(R£-R)  +  K^(V'f-^) 

+  K5(mf-m)2]  (66) 


Where  a  terminal  value  of  a  variable  xj^ 
is  unspecified,  the  corresponding  is 

taken  zero.  In  our  Illustrative  computa¬ 
tions,  we  have  chosen  the  problem  of  mini¬ 
mum  time  transfer,  that  is 


P  -  (67) 

We  note  that  the  function  H  ,  (64) , 

is  linear  in  i^,  a  feature  usually  asso¬ 
ciated  with  bang-bang  control.  As  recent¬ 
ly  pointed  out  by  Lawdenll,  however,  there 
is  a  possibility  that  the  collected  coef¬ 
ficient  of  q  in  H* 
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D  »  T2 


I  ^1  ^2 

—  sin  9  +  —  cos  9 
\  m  m 


V 

e 


may  vanish  over  a  finite  interval  of  time 
Such  arcs  satisfy  the  weak  form,  but  not 
the  strong  form,  of  the  Weierstrass  condl 
tion;  therefore  fall  in  the  gap  between 
necessary  and  sufficient  conditions  for  a 
minimum.  The  question  of  vjhether  an  arc 
D  =  0  may  or  may  not  be  minimizing  is 
currently  unresolved,  and  consequently  we 
have  no  a  priori  assurance  of  a  bang -bang 
throttle  characteristic. 


Some  Computational  Results 

The  three  successive  approximation 
schemes  described  earlier  were  mechanized 
for  digital  computation  of  planar  trans¬ 
fers  between  the  orbits  of  Earth  and  Mars, 
idealized  as  circular.  A  modified  Adams 
numerical  integration  scheme  was  used  with 
a  fixed  time  interval  of  two  days. 

Having  earlier  obtained  experience 
with  the  gradient  method  '.n  a  constant 
thrust  version  of  this  examplel®,  we  first 
performed  the  modification  to  incorporate 
the  throttle  variable  ri.  The  computer 
mechanization  employed  alternating  descent 
cycles  on  the  variables  9  and  rj.  After 
some  experimentation  to  select  penalty 
constants  K^,  a  family  of  transfers  for 
various  specified  terminal  values  of  the 
vehicle  mass  were  computed  (Figs.  1  and 
2)  .  Terminal  values  of  radius  and  veloci¬ 
ty  components  corresponding  to  the  Mars 
orbit  were  specified.  The  terminal  value 
of  the  heliocentric  angle  ij/  was  left 
open.  The  results  indicated  a  bang -bang 
throttle  characteristic,  the  transfers 
consisting  of  an  initial  full  throttle 
period,  a  coasting  period,  and  a  final 
full  throttle  period. 

One  such  transfer  was  adopted  as  a 
test  specimen  for  experimentation  with  the 
three  successive  approximation  techniques. 
The  initial  approximation  consisted  of 
circumferential  thrust  at  full  throttle  in 
all  three  cases.  Penalty  constants  were 
set  at  'intermediate'  values  for  the  first 
sixty  cycles,  then  increased  by  a  factor 
of  thirty.  It  was  intended  to  divorce  the 
effect  of  penalty  constant  manipulation, 
which  has  a  pronounced  effect  on  conver¬ 
gence,  from  the  characteristics  of  the 
successive  approximation  methods  in  this 
way.  This  was  not  very  successful  because 
of  strong  interaction  effects  —  the  pen¬ 
alty  constant  adjustment  technique  should 
be  'tailored'  to  the  method  and  to  the 
problem  under  attack. 


Successive  approximations  produced 
by  the  modified  gradient  method  are  shown 
in  Figs.  3A  and  36.  The  decrease  in  the 
function  p'  versus  number  of  descent 
cycles  is  illustrated  in  Fig.  3C. 

Corresponding  results  for  the  Min  H* 
scheme  are  shown  in  Figs.  4A,  4B,  and  4C. 
Increments  in  the  control  variables  9 
and  T]  were  generated  simultaneously  dur¬ 
ing  each  cycle.  Values  of  the  interpola¬ 
tion  parameter  C  corresponding  to  local 
minima  of  P  in  the  one -dimensional 
search  process  were  found  to  be  on  the 
order  of  .0003  to  .10,  the  smaller  values 
arising  in  conjunction  with  large  penalty 
constants . 

In  the  first  attempts  at  performing 
computations  with  the  integral -absolute 
value  metric,  a  difficulty  was  encoun¬ 
tered:  the  procedure  failed  to  continue 

to  produce  decreases  in  p'  after  six  or 
seven  cycles.  This  was  traced  to  the  fact 
that  the  search  for  a  variation  in  throt¬ 
tle  history  is  'quantized',  this  arising 
from  the  finite  number  of  numerical  inte¬ 
gration  steps.  Thus  the  'smallest  varia¬ 
tion'  in  throttle  history  admitted  by  the 
integration  procedure  is  lAp'  =  1  over 
one  integration  interval.  If  such  a  vari¬ 
ation  cannot  produce  a  decrease  in  P', 
the  one -dimensional  search  fails.  This 
difficulty  was  overcome  by  reducing  the 
integration  interval  to  one -tenth  of  that 
previously  employed,  i.e.  from  2  days  to 
.2  days.  The  increase  in  computing  time, 
by  a  factor  of  approximately  five,  indi¬ 
cates  a  need  for  employment  of  a  variable 
integration  step  feature  in  conjunction 
with  this  method. 

The  results  of  Figs.  5A,  5B,  and  5C 
were  obtained  using  alternate  cycles  of 
the  integral -absolute  value  scheme  on  q, 
and  the  Min  H*  scheme  on  6.  As  ob¬ 
served  earlier,  the  bang-bang  character 
of  the  throttle  history  is  preserved 
throughout  the  process,  this  being  a  fea¬ 
ture  of  the  method. 


Concluding  Remarks 

It  is  felt  that  the  numerical  results 
obtained  are  too  limited  to  provide  con¬ 
clusions  on  the  relative  merits  of  the 
three  methods.  The  differences  in  speed 
of  convergence  exhibited  in  Figs.  3,  4, 
and  5  are  insignificant  in  comparison  to 
improvements  attainable  by  modest  amounts 
of  experimentation  with  penalty  constant 
adjustment  procedures. 


Perhaps  the  most  significant  fact 
emerging  from  the  experiments  reported 
herein  is  that  several  successive  approxi¬ 
mation  techniques  can  be  successfully 
adapted  to  problems  featuring  inequality 
constraints  on  the  control  variables.  The 
three  methods  examined  are  workable  and 
possess  the  'hammer  and  tongs'  quality  so 
desirable  for  engineering  applications. 

On  the  other  hand,  it  is  clear  that 
continuing  research  in  the  class  of  suc¬ 
cessive  approximation  methods  is  likely  to 
be  fruitful,  perhaps  not  only  in  the  evo¬ 
lution  of  more  efficient  computational 
schemes,  but  also  in  contributing  to  the 
understanding  of  the  various  phenomena 
arising  in  variational  problems  of  flight 
performance,  control  and  guidance. 
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Symbols 

state  variables 

yj^  control  variables 

(■)  derivative  with  respect  to  time 

functional  representation  of  basic 
system  right  members  (Eq.  1) 

t  initial  time 

o 

tj,  terminal  time 

C)  initial  conditions 
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function  of  terminal  values  to  be 
minimized 

modified  function  including  penalty 
terms  (Eq.  3) 

positive  weighting  factors  for  pen¬ 
alty  terms 

state  variable  perturbations 
control  variable  perturbations 
function  defined  by  Eq.  8 
total  change  in  state  variable 
total  change  in  control  variable 
(0  <  <  1) ,  constant 

variable  for  nominal  trajectory 

I 

total  change  in  function  P 

system  variable  used  for  evalua¬ 
ting  AP' 

adjoint  variables  for  z  system 

function  defined  by  Eqs.  23  and  47 

total  change  in  y,  for  gradient 
direction 

positive  constant 

(0  >  bj^)  ,  search  parameter  for 
gradient  method 

lower  limit  on  control  variable  y|^ 
upper  limit  on  control  variable 

(0  <  C  <  1)  ,  searching  parameter 
for  Min  H*  method 

"k 

control  function  for  Min  H 

search  parameter  for  Integral - 
absolute  value  metric  method 

generic  functions 

radial  velocity 

circumferential  velocity 

radial  distance 

thrust 


gravitational  constant 
m  mass 

9  thrust  direction  angle 

f  heliocentric  angle 

exhaust  velocity 

q  throttle  control  variable 

Aq  increment  in  throttle  control 

variable 

(  )^  final  values 

k 

D  coefficient  of  q  in  H 
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FIG.  3A  SUCCESSIVE  APPROXIMATIONS  TO  PLANAR  EARTH -MARS 
ORBIT  TRANSFER  VIA  MODIFIED  GRADIENT  METHOD 
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FIG.  3B  SUCCESSIVE  APPROXIMATIONS  TO  PLANAR  EARTH-MARS 
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FIG.  4A  SUCCESSIVE  APPROXIMATION  TO  PLANAR  EARTH -MARS 
ORBIT  TRANSFER  VIA  MIN  H*  METHOD 
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FIG.  5A  SUCCESSIVE  APPROXIMATIONS  TO  PLANAR  EARTH -MARS 

ORBIT  TRANSFER  VIA  INTEGRAL  -  ABSOLUTE  VALUE  METHOD 
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FIG.5B  SUCCESSIVE  APPROXIMATIONS  TO  PLANAR  EARTH- MARS 
i50r  ORBIT  TRANSFER  VIA  INTEGRAL-ABSOLUTE  VALUE  METHOD 
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